﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: ObjectiveFunctionBase.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2013-01-24 4:34 PM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using System.Linq;

namespace BOL.Algorithms.Optimization
{
    public abstract class ObjectiveFunctionBase : IObjectiveFunction
    {
        #region Private variable

        private double[] _current;
        private List<Predicate<double>>[] _constraints;
        private const double Epsilon = 1e-8;

        #endregion

        #region Public properties

        /// <summary>Gets the current parameters of objective function.</summary>
        public double[] Current
        {
            get { return _current; } 
            set
            {
                if (value == null)
                    throw new ArgumentNullException("value");

                if (value.Length != NumberOfParameters)
                    throw new ArgumentException();

                _current = value;
            }
        }

        /// <summary>Gets the number of parameters of objective function.</summary>
        public int NumberOfParameters { get { return _current.Length; } }


        public List<Predicate<double>>[] Constraints
        {
            get { return _constraints; }
            set
            {
                if (value == null)
                    throw new ArgumentNullException("value");

                _constraints = value;
            }
        }

        #endregion

        #region Constructor

        protected ObjectiveFunctionBase(double[] initial)
        {
            if (initial == null)
                throw new ArgumentNullException("initial");

            _current = initial;
            _constraints = initial.Select(x => new List<Predicate<double>>()).ToArray();
        }

        protected ObjectiveFunctionBase(int numberOfParameters)
            : this(new double[numberOfParameters]) { }

        #endregion

        #region Public methods

        /// <summary>Gets the value of objective function at the current position.</summary>
        public double Value()
        {
            return Value(_current);
        }

        /// <summary>Gets the value of objective function at a specific position.</summary>
        public abstract double Value(double[] parameters);

        /// <summary>Gets the gradient of objective function at the current position.</summary>
        public double[] Gradient()
        {
            return Gradient(_current);
        }

        /// <summary>Gets the gradient of objective function at a specific position.</summary>
        public virtual double[] Gradient(double[] parameters)
        {
            if (parameters == null)
                throw new ArgumentNullException("parameters");

            if (parameters.Length != NumberOfParameters)
                throw new ArgumentException();

            var gradient = new double[NumberOfParameters];

            for (var i = 0; i < NumberOfParameters; i++)
            {
                parameters[i] -= Epsilon;
                var fm = Value(parameters);
                parameters[i] += 2 * Epsilon;
                var fp = Value(parameters);
                gradient[i] = .5 * (fp - fm) / Epsilon;
                parameters[i] -= Epsilon;
            }

            return gradient;
        }

        /// <summary>Gets the hessian of objective function at the current position.</summary>
        public double[][] Hessian()
        {
            return Hessian(_current);
        }

        /// <summary>Gets the hessian of objective function at a specific position.</summary>
        /// <seealso cref="http://en.wikipedia.org/wiki/BFGS#Algorithm"/>
        public virtual double[][] Hessian(double[] parameters)
        {
            if (parameters == null)
                throw new ArgumentNullException("parameters");

            if (parameters.Length != NumberOfParameters)
                throw new ArgumentException();

            var hessian = new double[NumberOfParameters][];

            for (var i = 0; i < NumberOfParameters; i++)
            {
                hessian[i] = new double[NumberOfParameters];
                parameters[i] -= Epsilon;
                var fm = Gradient(parameters);
                parameters[i] += 2 * Epsilon;
                var fp = Gradient(parameters);
                for (var j = 0; j < NumberOfParameters; j++)
                    hessian[i][j] = .5 * (fp[j] - fm[j]) / Epsilon;
                parameters[i] -= Epsilon;
            }

            return hessian;
        }

        #endregion
    }
}
